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Studies of past sea-level markers are commonly used to unveil the tectonic history and seismic behavior of 
subduction zones. We present new evidence on vertical motions of the Hellenic subduction zone as resulting 
from a suite of Late Pleistocene - Holocene shorelines in western Crete (Greece). Shoreline ages obtained by 
AMS radiocarbon dating of seashells, together with the reappraisal of shoreline ages from previous works, 
testify a long-term uplift rate of 2.5-2.7 mm/y. This average value, however, includes periods in which the 
vertical motions vary significantly: 2.6-3.2 mm/y subsidence rate from 42 ka to 23 ka, followed by 
—7.7 mm/y sustained uplift rate from 23 ka to present. The last ~5 ky shows a relatively slower uplift rate of 
3.0-3.3 mm/y, yet slightly higher than the long-term average. A preliminary tectonic model attempts at 
explaining these up and down motions by across-strike partitioning of fault activity in the subduction zone. 

Studies of past sea-level markers are commonly used to unveil the tectonic history and infer the seismic 
behavior of subduction zones 1 " 3 . These studies usually provide estimates of upper-plate long-term uplift 
trends often with a large variability at the 10 5 -year timescale 4 ' 5 . Explanations proposed for this variability 
include the roughness of the subducting plate and down-dip variations of the locked zone 6 . Vertical motion also 
manifests sign changes at the timescale of 10 2 -10 3 years due to the subduction earthquake cycle 7 , and at few-years 
distance on occasion of distinct earthquakes 8 . 

In the Mediterranean Sea, the Hellenic subduction is the one with a supposed capability of generating very large 
earthquakes and tsunamis. This awareness is mainly brought forward by the historical notion of the 365 AD 
event 9 " 11 . In this subduction zone, the island of Crete is one of few accessible places that may provide tectonic 
information on the overriding plate (Figure 1). In the last 5 My, the area has been uplifted by —1000 m 12 . In 
recent times, uplift is testified by raised shorelines and marine terraces in southwestern Crete reported in studies 
as old as that of Spratt 13 and recent ones, e.g. Peterek et al. 14 (and references therein). Nearby Paleochora, several 
shorelines are recognizable at various elevations; some of them are of Late Pleistocene age and correlated with 
eustatic peaks (Supplementary Table SI online). Using these data, the average uplift rate for the last 40-50 ky is 
estimated at ~ 1.5 mm/y by Wegmann 15 and ~2 mm/y by Shaw et al. 16 . Another estimate at a location 33 km east 
of Paleochora is of 1-1.5 mm/y 17 . 

Besides this net uplift trend, a first indication of both upward and downward movements in Crete is suggested 
by ten shorelines progressively younger from bottom up and dated between 4800 y BP (the lowermost) and 
1550 y BP (the highest) that testify ten small (10-20 cm) lowering events, for a total of about 3 m 18 . These 
shorelines then became exposed by the up to ~9-m uplift ascribed to the 365 AD earthquake by Pirazzoli et al. 19 , 
Shaw et al. 16 , Stiros and Papageorgiou 20 , and Stiros 21 . 

In search for additional constraints on the uplift history of the Hellenic subduction, we surveyed the Paleochora 
shorelines, constrained their age through AMS radiocarbon dating, and then combined these data with paleo sea- 
level reconstructions for the past —50 ky. 

With respect to previous estimates, our results suggest a faster net uplift rate and 10 4 -year-long periods of 
alternating uplift and subsidence, during which the vertical movements are very different from the net average. 
Although a thorough tectonic explanation of these vertical movements is beyond the scope of this work, this 
behavior may be explained by across-strike partitioning of tectonic activity (e.g. several subduction earthquake 
cycles) in the subduction zone. Knowledge of these variations can contribute to better understanding the long- 
term behavior of the Hellenic subduction zone, and help constrain its earthquake productivity. 

Results 

Raised shorelines. The coast of Crete between Krios and Paleochora (~8 km) preserves a rich set of rather 
continuous raised shorelines (Figure 2 and Figure 3). Apart from the stretch between Krios and Plakaki, alluvial 
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Figure 1 | Tectonic sketch map of the Hellenic subduction zone. The black square indicates location of the studied coastline nearby Paleochora. 
Topographic and bathymetric relief are obtained using SRTM V2 (http://www2.jpl.nasa.gov/srtm/) and GEBCO (http://www.gebco.net/) data, 
respectively. 

fans and deep incisions alters the lateral continuity of shorelines only meters (Figure 3; Supplementary Table S2 online), of Late 
locally. We here illustrate six shoreline remnants, named SI through Pleistocene - Holocene age (Table 1; Supplementary Table S3, 
S6 from bottom up, irregularly spaced at elevations between 4 and 75 Figure SI, Figure S2 online). 




Figure 2 | Pictures of the surveyed shorelines. Kalamia site notches and their elevations (see Figure 3 for location): (a) Panorama showing three raised 
shorelines found in the West Kalamia site; (b) A band of lithophagid boreholes marking shoreline S2 in East Kalamia site; (c) S5 shoreline in the West 
Kalamia site; notice the cave in the background and the sandy-gravel deposit on the foreground; (d) S6 shoreline in the West Kalamia site; notice the cave 
on the left and the notch carved in the bedrock, filled with gravel. 
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Figure 3 | Paleochora shoreline suite. Sketch profile of the whole Kalamia site (left) and synoptic view (right) of all the raised shorelines and their ages in 
the Paleochora broader area; major and minor shoreline features are indicated by thick and thin segments, respectively; vertical extent of sedimentary 
bodies is indicated by double arrows. Other symbols: @ = fossil rinding; dotted line = shoreline correlation. 



SI is very continuous and well exposed and it is actually composed 
by a suite of several features of different ages. The highest in the suite, 
denoted SI -high, consists of a notch with algal concretions at 8.34 m 
and a surf notch reaching 10 m. At places, remnants of a pebbly 
beachrock adjoin the main notch, and one such pebbly beach deposit 
at 8.8-11.8 m of elevation includes a layer rich in bivalve shells on 
apparent lateral continuity with an abrasion platform. Several other 
minor notches are present below the main notch at elevations of 4- 
8 m. We denote the lowermost notch in this suite as SI -low. AMS ages 
allow correlation of the entire SI suite with that described by Pirazzoli 
et al. 18 in other parts of Crete (Supplementary Table S2 online). 

S3 is located at 16.5 m and is also very well exposed. It can be 
followed almost continuously for hundreds of meters and is charac- 
terized by a well developed notch and lithophagid boreholes, with 
adjoining sea arches, caves, and small abrasion platform remnants. 
Several minor notches can also be seen below the main feature. S3 
also features a — 7-m-thick beach deposit that yielded a bivalve shell 
of mid-Holocene age. 

S2, S4, and S5 are located at elevations of 14, 34, and 55 m, respect- 
ively. The age of S2, whose remnants are essentially lithophagid 
boreholes and a notch, is constrained by a Lithophaga sp. shell dated 
21-22 ka. The S4 age of 40-48 ka is consistently provided by four 



samples. Two of them were found included at the top and bottom of a 
1.6-m-thick algal-reef remnant in the Paleochora peninsula, which is 
geometrically correlated with a gravel deposit at 34 m elevation on 
the coastal slope that yielded other two bivalves of about the same 
age. Conversely, the three samples collected in the gravel associated 
with S5 are of different ages spanning a considerable time interval of 
22-44 ka. The oldest sample was found in loose sediment and can 
thus be suspected of having been misplaced. S6 is poorly preserved; it 
only has a small exposure of a notch at 75 m with sterile gravel and 
could not be dated. 

The analyzed series of shoreline remnants is entirely exposed only 
at Kalamia. Similar shoreline remnants can also be seen in other parts 
of Crete but they are not considered in this study. 

Average long-term uplift. Sea-level variations for the period of 
interest (—50 ky), which includes the last three Marine Isotope 
Stages (MIS), are well represented in the eustatic curve by Wael- 
broeck et al. 22 (Supplementary Figure S3 online). However, consider- 
ing the scatter among curves in different periods 23 , resorting to 
specific literature data for each single time interval is necessary. 
Most authors agree on a — 120 m sea level during the MIS2 
(—20 ka). A recent review for MIS3 24 indicates a sea level of 
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Table 1 | Sample dating results 



Code 


Shoreline and Locality 


Description 


dl3C 


Conv. AMS Age (y BP) 


2a Cal. Age (y BP) 


WKB1B 


SI; WKalamia 


Vermetid shell 


-0.5 


2190 ± 0.040 


1519-1933 


WKB1 A 


SI; WKalamia 


Spondylus sp. 1 


-0.9 


2260 ± 0.040 


1591-2029 


WPAB41 


SI; P . Ammos 


Bivalve shell 


-8.4 


2348 ± 0.033 


1 707-2 1 25 


WKB4 


SI ; W Kalamia 


Lithophaga sp. 
Bivalve shell 2 


2.5 


2710 ± 0.040 


2142-2652 


KB26A 


S3; Krios 


1.6 


5162 ± 0.039 


5286-5639 6 


EKB23 


S2; E Kalamia 


Lithophaga sp. 


-7.6 


1 8400 ± 0.080 


21 040-21 880 


WKB5 


S5; W Kalamia 


Bivalve fragment 


-3.1 


19480 ±0.090 


22320-23030; 23070-23290 6 


WKB13 


S5; W Kalamia 


Bivalve fragment 


0.1 


26930 ±0.160 


30810-31300 


CA43 


S4; Paleochora 


Bivalve fragment 


-8.9 


36993 ± 0.642 


40340-42530 


WKB10 


S4; W Kalamia 


Bivalve fragments 3 


-1.5 


37530 ±0.330 


41380-42450 


WKB8 


S4; W Kalamia 


Glycymeris sp. 2 


0.2 


38250 ±0.350 


41840-42950 


WKB6 


S5; W Kalamia 


Bivalve fragment 4 5 


-0.5 


39190 ±0.380 


42370-43770 


CA42 


S4; Paleochora 


Bivalve shell 


-3.5 


42402 ±1.102 


43490-47430 



Notes: Results are rounded to the nearest 1 0 y for samples with standard deviation in the radiocarbon age greater than 50 y. 
Wg + Cal + Qz; 
2 Arg + Cal < 5%; 

3 Four fragments: 1 ) Arg; 2) and 3) Arg + Cal; 4) no XRD; 
4 Arg; 

5 Loose fragment at the base of a conglomerate deposit; 
6 Full range of multiple calibration peaks. 



— 60 m in the early part (—50-60 ka) and of — 80 m in the late part 
(—30-50 ka). For the Holocene in Crete, two models exist: 1) same 
sea level elevation as of today 25 , and 2) sea level 1.25-1.5 m lower 
than present 26 at —6 ka. 

With respect to these past sea-level elevations, the Paleochora 
shorelines collectively indicate a net long-term (—45 ky) uplift rate 
for western Crete of 2.5-2.7 mm/y (Figure 4). 

Ups and downs. Are some of the shoreline points in Figure 4 full- 
fledged indicators of important deviations from the average long- 
term uplift trend? 

Consider first that the age of S4 can be placed in the late MIS3, 
when sea level was at —80 m, as constrained by four samples with 
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Figure 4 | Vertical displacement vs. shoreline ages. Yellow dashed line 
represents the long-term average uplift trend; blue dashed line depicts 
main up-and-down trend. Notice that Kell4C and Shal4C ages are not 
calibrated. The star symbol marks a sample age discarded by Wegmann 15 
because of contamination. Further details in Supplementary Tables SI and 
S3 online. 



consistent ages (WKB8, WKB10, CA43, and CA42). The oldest age of 
S5, instead, cannot be accurately traced because of the uncertain 
placement of WKB6 and sea-level oscillations as large as 20 m during 
MIS3, but one can be confident in considering that WKB13 and 
WKB5 imply that S5 was active in the late MIS3 and early MIS2 with 
sea level lowering from —80 m to — 120 m. 

Attaining the S4 to S5 vertical separation thus requires a net sub- 
sidence rate of 2.6-3.2 mm/y in the period from —42 to 23 ky ago. A 
period of sustained uplift of —7.7 mm/y should have then followed 
the S5 abandonment (—23 ky ago) as suggested by the formation of 
S2, with sea level at —120 m, and S3 and SI -low with sea level at 
about the same elevation as today. Episodic downward movements 
led to the formation of the other several shorelines of the SI suite and 
SI -high itself, which were then all displaced by the uplift event in 365 
AD 19,21 . For what concerns S3, finding of KB26A partly fills the gap of 
lower-mid Holocene data pointed out by Kelletat 27 and provides the 
first evidence for a Holocene shoreline above SI -high. S3 should have 
been abandoned in a period between the age of KB26A and the age of 
the oldest shoreline of the SI suite (Figure 3). From our determina- 
tions (Table 1) Sl-lowhas an age of 2142-2652 y BP, thus suggesting 
the S3 abandonment to have occurred within 2142-5639 y BP. This 
time range, however, could likely be much shorter if considering that 
SI -low is younger than shoreline VIII 19 , found at 4 m in Phalasarna 
(Crete western coast), whose age is 4410-5070 y BP 19 . The net late- 
Holocene uplift rate in Crete is slightly different depending on which 
of two sea-level models is adopted for the time of S3 abandonment. 
The uplift rate would be 3 mm/y according to sea level as in the work 
by Pirazzoli 25 ; whereas it would be 3.3 mm/y according to sea level as 
in the work by Lambeck and Purcell 26 . 

Preliminary tectonic model. Mature accretionary convergent 
margins are known to develop distinct kinematic domains across- 
strike which can also host splay faults. Along-dip partitioning of the 
slab itself are also documented in subduction zones and have been 
explained by various driving mechanisms 28 . In the Hellenic 
subduction, the best known uplift event in Crete, i.e. the event 
associated with the 365 AD earthquake, has either been related to 
seismic slip on the subduction interface 29,30 or on a shallower crustal 
fault embedded into the accretionary wedge 16,21 . Building on this 
knowledge and as a preliminary attempt to explain the 
reconstructed vertical motions, we envisage that 10 4 -year timescale 
rate variations can be explained by fault activity partitioning across 
the subduction zone. To this end, we subdivide the subduction zone 
into two separate sections, front and rear, across the main slab dip 
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Figure 5 | Fault dislocation model. Steady-state (t -> infinity) surface 
displacement pattern (Uz) along profile shown in Figure 1 due to faulting 
in a viscoelastic-gravitational medium 43 . We hypothesize four equal-size 
fault ruptures (roughly moment magnitude 8) at the front and rear of the 
slab main dip change (vertical exaggeration —2.5 X). Front section: Al, 
interface (solid blue); A2, splay (dashed blue). Rear section: Bl, interface 
(solid red); B2, splay (dashed red). More modeling details are shown in 
Supplementary Tables S5, S6, and S7 online. Notice that the relative fault 
rupture positions and the fault dip angles determine vertical movements of 
opposite sign at the observation point (Paleochora). 

change and set up a series of four distinct dislocation models which 
include fault slip on the subduction interface and on prospective 
splay faults (Figure 5). 

Based on these fault- dislocation models, we suggest that periods of 
enhanced uplift rates at Paleochora occur when activity in the sub- 
duction rear section dominates (i.e. when the fault activity mainly 
occurs beneath Crete, as in the case of the 365 AD earthquake). When 
considering splay ruptures through the upper plate the uplift trend is 
further emphasized. Because of their high dip angle these faults 
mainly produce uplift and limited subsidence. Conversely, stability 
or lowering at Paleochora occurs when the fault activity primarily 
involves the subduction front section (i.e. south of Crete) including 
shallow splay faults within the accretionary wedge. Persistence of this 
activity can explain both the lowering episodes recorded by the SI 
suite 18 and the subsidence between S4 and S5. Because of their low 
dip angle, these faults produce significant subsidence in Crete and 
uplift offshore south of Crete. 

The alternating phases of uplift and subsidence recorded by the 
raised shorelines on the coast of Crete may thus reflect backward and 
forward shifting of tectonic activity across the Hellenic subduction 
system as also observed at a different timescale in other subduction 
zones (e.g. Solomon Islands 8 ). 

Discussion 

Major uncertainties in our reconstruction of the Paleochora shore- 
line suite are represented by the measure of their elevation and spatial 
correlation, and the sample age determination and sampling site 
conditions. Shoreline elevations are measured by a geodimeter whose 
accuracy, together with the small tidal range in Crete, ensure an 
overall elevation uncertainty of well developed notches to remain 
within less than ±10 cm. Uncertainty can be occasionally higher 
for gravel deposits or other coarse shoreline features. The short hori- 
zontal distance that separates most remnants allows us considering 
the lateral correlation based on geomorphic criteria to be relatively 
robust. We cannot exclude that the algal reef in the Paleochora pen- 
insula can be correlated with S5 but this alternative would not affect 
our interpretation of the results. X-ray diffraction analysis shows that 
radiocarbon ages are not affected by calcite recrystallization. The 
small percentage of calcite (<5%) in KB26A, WKB8, and WKB10 



(Supplementary Figure S2, Table S3 online) may indicate secondary 
recrystallization but not necessarily alteration of their radiocarbon 
age. As for KB26A, the effect of contamination would result in reju- 
venating its age of at maximum 200 y. Accordingly, the older bound 
of its calibrated age may be extended to —5880 y BP. Considering 
that real ages of Lithophaga sp. shells can be up to 2 ky younger than 
their radiocarbon age 31 , the younger bound of EKB23 calibrated age 
could be up to ~ 19 ka. The site conditions (Supplementary Table S2 
online) indicate a strong bond between samples and shoreline in all 
cases except for WKB6, suggesting that its age is not totally reliable to 
ascertain the age of S5. 

The accuracy of uplift rate estimates is affected by all the above 
uncertainties along with uncertainties about past sea-level eleva- 
tions 23 . Our net long-term (—45 ky) uplift rate estimate of 2.5- 
2.7 mm/y (Figure 4) is higher than the uplift rate (~2 mm/y) 
obtained by Shaw et al. 16 which consider an elevation of the 
Paleochora terrace of 24 m instead of the inner edge elevation of 
34 m resulting from our survey (Supplementary Table S2 online). 
The average uplift rate proposed by Wegmann 15 , instead, is signifi- 
cantly lower (—1.5 mm/y) than our estimate. We note that this value 
is mainly based on the correlation of a shoreline at 9 m elevation 
(likely SI -high in our work) with MIS3. However, these differences 
are relatively small compared with the uncertainty in sea-level eleva- 
tion at the time of shoreline formation or abandonment and mainly 
depend on interpretations of individual shorelines. Nonetheless, the 
reconstructed pattern of ups and downs clearly deviates from the 
average trend (Figure 4) supporting the idea that vertical displace- 
ment in Crete may strongly fluctuate over a 10 4 -year timescale. For 
completeness of information, we note that some other shoreline ages 
nearby Paleochora (Supplementary Table SI online) are difficult to 
reconcile with our reconstruction. For example, Kelletat and 
Zimmermann 32 provided an uncalibrated age of 5.735 ± 90 ky BP 
for an algal rim at 3.3 m, at Kalamia, 13 m below S3; Kelletat et al. 33 
determined an ESR age of 113-114 ka for two samples collected at 
2.5 m elevation; and Wegmann 15 determined a radiocarbon age of 
41-46 ka for the shoreline at 9 m elevation. However, Wegmann 15 
also provided 36 CI exposure ages for some raised shorelines around 
Paleochora that, when combined with the eustatic curve, fit quite well 
with our reconstruction (Figure 4). 

Because of the dominating role of convergence in the Hellenic 
subduction, our preliminary tectonic model only focuses on fault 
activity in the subduction zone and neglects possible contributions 
from other processes, such as the crustal extension in the upper plate 
(see for example modeled geodetic rates from Reilinger et al. 34 ). Nor 
have we investigated the role of the various driving mechanisms of 
fault activity in the subduction system (see Kopp 28 for a review). 
Nonetheless, at this stage of the analysis we propose a simple model 
that together with our findings shows that tectonic rates may vary 
depending on both timescale and location of the observation point 
with respect to the subduction architecture. As we progressively 
unveil vertical tectonic rate fluctuations at various timescale in sub- 
duction zones, the analysis of raised shorelines proves to be an effec- 
tive tool to improve our understanding of long-term processes which 
complement other observables such as decadal instrumental mea- 
surements. In this perspective, detailed shoreline analyses can shed 
light on maximum deviations of tectonic rates with respect to long- 
term averages. Since earthquake productivity estimates can be derived 
from tectonic rates (e.g. Geist and Parsons 35 ), a compelling implica- 
tion in active subduction zones is that earthquake rates may consid- 
erably vary as a function of location and time interval considered. 

Methods 

Shoreline identification. Raised shorelines are identified by coastal notches, 
lithophagid boreholes, sea caves, coastal terraces, abrasion platforms, and sand-to- 
gravel beach deposits. Their present elevations are measured through geodimeter 
land surveying. Correlation of individual shoreline remnants is initially based on 
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geomorphic criteria and then confirmed/refined through AMS radiocarbon ages 
where applicable. 

Shoreline ages. The AMS radiocarbon technique is used for dating seashells or 
encrusting organisms sampled directly from the shoreline feature. X-ray diffraction 
analysis (XRD) of the inner parts of dated samples is used to verifying the 
preservation of the original aragonite phase of shells 36 . Pretreatment of samples 
includes mechanical and chemical removal of the outer portion in order to limit the 
effects of external contamination. Calibration of radiocarbon ages is carried out using 
INTCAL09 Marine Calibration dataset 37 with the Calib 6.0 software (© 1986-2011: 
Minze Stuiver and Paula Reimer) and, following recommendations by Reimer and 
McCormac 38 , by taking into account the local reservoir offset from the modeled world 
ocean (AR) by applying a value of AR = 46 ± 77 y as the Mediterranean average 
(Marine Reservoir Correction Database, http://calib.qub.ac.uk/marine/, last accessed 
March 3rd, 2014; using data from various sources 38 " 41 , including a pers. comm. to the 
authors of the Database by Taviani and Correggiari; Supplementary Table S4 online). 
Laboratory error multiplier equal to 1 has been applied to all samples. When 
discussing radiocarbon ages taken from the literature we refer to values that are 
calibrated or re-calibrated in this same way. 

Shoreline displacement rates. Radiocarbon ages are used for the chronological 
correlation of shorelines with the eustatic curve and MISs (Supplementary Figure S3 
online). We assume that a shoreline or terrace may form during any time when sea 
level and tectonic variations were about equal 42 , then obtain the net shoreline 
displacement by subtracting the eustatic elevation at the time of formation/ 
abandonment from the measured shoreline elevation. Vertical displacement of 
shorelines is defined as AU = H — AE, where: A U is the calculated displacement, H is 
the measured shoreline elevation above present sea level, and AEis the difference 
between present mean sea level and sea-level elevation (eustatic) at the estimated time 
of shoreline formation. Vertical displacement rates (either uplift or subsidence) are 
obtained by averaging the ratio between vertical separation of shorelines and the 
intervening time in a period of interest, such as UR = (AU n — AU n -i)/T n — T n - 1 
where: UR is the calculated displacement rate in the time interval T n — T n -i\AU n and 
T n are the calculated displacement and the age of shoreline n, respectively. 

Preliminary tectonic modeling. To explore the vertical component of the surface 
displacement patterns along a profile normal to subduction strike and through 
Paleochora (Figure 5) we make use of viscoelastic-gravitational fault-dislocation 
modeling 43 . All faults are purely reverse, and strike orthogonally to profile in Figure 1 
with dip toward NE according to the local slab geometry 44 ; they have the same size and 
amount of slip as to represent a magnitude- 8 earthquake compatible with appropriate 
scaling laws for slab interfaces 45 and crustal faults 46 . Modeling setup details are shown 
in Supplementary Tables S5, S6, S7 online. Supplementary Figures S4 and S5 online 
also show profiles of both coseismic and steady-state displacement for the front and 
rear sections of the subduction zone. 
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